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ABSTRACT 

We propose a new method to calculate expectation values of a delta function of 
the Hamiltonian, (\1' | 5{H — E) \ \E'). Since the delta function can be replaced with 
a Gaussian function, we evaluate 

(VP I ^e-'^(^-^)' I ^!) 

with large j3 adopting the Suzuki- Trotter decomposition. Errors of the approximate 
calculations with the finite Trotter number Nt are estimated to be 0{1/N^) for 
the i^th-order decomposition. The distinct advantage of this method is that the 
convergence is guaranteed even when the state | \E') contains the eigenstates whose 
energies spread over the wide range. 

In this paper we give a full description of our method within the quantum mechan- 
ical physics and present the numerical results for the harmonic oscillator problems 
in one- and three-dimensional space. 
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1 Introduction 



In the quantum physics some dynamical quantities in the matter[|I|, |^] can be expressed by 

a{E) = Im{(* I A''-. ^- A | *)} = 7r(* | - E)A \ (1) 

H — E — ie 

Examples of them are the density of states and the forward scattering amplitude. There are many 
attempts to calculate the delta function by the polynomial expansions such as the Chebyshev 
polynomial expansion Q. In these expansions it is assumed that the eigenvalues are situated inside 
the fixed ranges. Whether one can obtain a good convergence or not with the expansion, therefore, 
entirely depends on the state | ^E"). 

In this paper we investigate a new method to numerically calculate the expectation value of the 
delta function. Our starting point is that the delta function can be represented by the Gaussian 
function 

I A^5{H - E)A I = lim I A\ (le~'^^"-'^'>' A \ (2) 

It is justified in the evaluation to use y^f3/n ■ e~^^^~^^'^ with large but finite (3 instead of the 
delta function because can be interpreted as the resolution in the actual observations. We 
then employ the Suzuki- Trotter decomposition to calculate the expectation value of the Gaussian 
function. This decomposition is quite stable as is well known in the studies of the quantum Monte 
Carlo method ||. Extensive work has also been done to calculate e by this compositionH. 

In the next section we explain our method in detail, applying it to the one-dimensional quan- 
tum mechanics. In Sec. 3 the higher-order decompositions Q are described and Sec. 4 is devoted to 
estimating errors in the method. We then present numerical results for the one-dimensional har- 
monic oscillator in Sec. 5, where the stability of the method is also demonstrated. Sec. 6 contains 
applications to the three-dimensional problems and the final section is for the summary. 

2 Methods 

Our study here will be limited to the quantum mechanics, i.e. to the one particle problems. Here 
we describe the method in the one-dimensional case. The Hamiltonian is given by 

H--^ + Vi.), (3) 

where we adopt units of h — 1 and m = 1/2. From Eq.(||) it follows 



^ The Gaussian function has been used for the filter function, where the polynomial expansion and other methods 
are employed to calculate itM. 
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Next we employ the discrete space representation, 

Xi — {i ~ l)-^ + Xmin — Ij ' ' ' : L) , A — iXmax ~ Xmin) / L. (5) 

The wavefunction (pix) — {x \ <j}) is replaced by 4>{xi)^ which is denoted as (pi hereafter. Then the 
differentials become 

^-^(0.+i+0.-i-20.), (6) 
d^Mx) 1 

^ ^(<^»+2 + </'»-2 - - Hr^l + (7) 

^ d^jx) , cfV{x)<ly{x) 

-^{[V{x,+i) + + [ V{xi) + - 4y(a;.)'/'J, (») 

respectively, and we obtain a matrix _ff2 to represent [H — i?)^. In order to calculate e~^^^ using 
the Suzuki- Trotter decomposition, we divide H2 into four matrices, i?2"^ (^^ — 1; 2, 3, 4), which are 
formed by aligning the 4x4 matrices along the diagonal line (See Fig. 1), 

In the second-order decomposition we approximate e~^^^ by 

exp(-^i?(^))exp(-^i?f' )ea:p(-^i?f))exp(-^i7(^))]^' (10) 
with the finite Trotter number Nf. The error in this approximation is 0{1/Nf), as is well known. 



3 Higher-order decompositions 

In Ref. the higher-order decompositions have been discussed. This section is to briefly rcfcr, for 
later use in our calculations, to the fourth- and the sixth-order decompositions together with the 
second-order one. For simplicity we describe here only the case where the operator C is the sum 
of the two operators A and B. 

The iiTth-order formula is defined by 

g.C ^ ^ [SK{x/Nt)f* + 0(a;^+ViVf ). (11) 

For the second-order decomposition {K — 2) we have 

^2(x) =e*^e"^et^. (12) 
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Using 5*2 the fourth-order decomposition {K — 4) is given by 

Siix) = S2{p2xfS2{{l - ^P2)x)S2{p2xf, (13) 

with p2 = {4- 4^/^)-^ = 0.4144907717943757, while the sixth-order one (A' = 6) is represented as 

Seix) = S2{pix)'^ S2{P2X)'^ S2{P3X)'^ S2{P4X)'^ S2{P3X)'^ S2{P2X)'^ S2{pix)'^ (14) 

with Pi = 0.3922568052387732, p2 = 0.1177866066796810, 
P3 ^ -0.5888399920894384 and pi = 0.6575931603419684. 

Since we have aheady obtained 5*2 for our problem in the previous section, it is straightforward 
to apply these higher-order decompositions to our calculations. 

4 Errors in the Suzuki- Trotter decomposition 

In this section we roughly estimate errors in our method using eigenstates | ■0^) and eigenvalues Xi 
for the hermitian operator C, namely 

C I ^P,) =1 V.)A„ (15) 

with {^i I ipj) ^ 6ij. 

For the operator S'_r-((5) in the Xth-order formula, where S denotes some small c- number, we 
have 

SK{S) = e'^ + Ek, 

where Ek is 0{S^^^). Then we apply the decomposition for e'^^*'^ with a fixed SNt, 

Nt~l 

= e^'^^+ J] e'=^'^^ife(^'-^-i)*'5 + 0((5Wi)). (16) 

k=0 

Note that the S should be 0{1/Nt) here. In order to go on with the estimation we neglect the oper- 
ators of 0((52(-^+i)) in Eq.(|l|), although this is not a rigorous procedure. With this approximation 
we obtain 

Nt~l 

= I e^'*^- -f ef^'^^ Ekc^^'-"-'^'^' I V.) 

fe=0 

Nt-l 

= + e'=^^^e(^'-'^-i)*^'(V, | Ek \ i^i). (17) 

k=0 
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For j — i it reads 

(7^, I [Sk{S)]''' I ^P^) 0, e^'*^'[l + I Ek I 4'^)Nte-'^^]. (18) 
When j ^ i, on the other hand, we have 

k=Q 

NtSX, _ NtSXj 

= {^^\Ek\^.) ^sx._^sx, ■ (19) 
Eq.(|l9|) suggests that the error is - | i^^ I ^p^)Nte'^^^~^'^^^^ - 5^ if 6{Xj - A,) « 0, while for 



the case > A, it is estimated by (V'j | Ek \ il)i)e^^'-^^^^^ - 5^'+^ 



5 Results in one dimension 

Now we present our results for the one-dimensional harmonic oscillator whose potential is 

F(a:)=Ax^ (20) 

where we fix A = 1/4 so that the energy eigenvalues are (fc— 1/2) (fc = 1, 2, 3, • • •) in the continuum 
limit. In our calculations with the discrete space representation using Eqs. (|^) ~ (|8|) we scale the 
matrix H2 by A^, which demands that the energy should be A^(fc — 1/2) instead of (fc — 1/2). 
Throughout this section we set Xmin — —25, x^ax — 25 and L = 500 in Eq. (|^). 

Let I (f)'-'^')) be the eigenstate of the Hamiltonian with the eigenvalue Ek, for which the ordering 
Ek < Ek+i is assumed. First, following to Eq.(p^ and employing the fourth-order decomposition, 
we estimate the errors of I e-'^^^--^)" I 0(20)^ ^ii-i^ p ^ iQQ foj. several states | In order 

to estimate errors owing to the decompositions only, we compare our numerical results with those 
obtained by the exact diagonalization of the matrix H2. In Fig. 2 we plot the results as a function 
of the Trotter number Nt, which clearly indicate that the errors decrease as 0{N^^). Fig. 3 is to 
compare the errors estimated with Nt — 1000 and those with Nt = 2000 for each k up to 40. Here 
we see that we can make the errors quite small by increasing Nt- Note that the largest error is 
found at fc = 20, which is the case described by the Eq. (^^. 

Next we calculate 

/(i;)^(0|e-«^-^)^|^), (21) 
using the wavefunction (pi which is parametrized by 

</>,: = </>r^Ci + ... + 4>f^CL, Cl + ... + = 1, (22) 
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with the discretized wavefunctions 0^'' =| (l)^^\xi)). Here we use two wavefunctions as (pi. One of 
them, which we will call the wavefunction (A), is given by 

Cfe = ^= for k<20, Ck = for k > 21. 
V2O ~ 



For the wavefunction (B), which is useful to study the effects due to the highly excited states, we 
set 

Ck - + — ) for k < 20, 

^V20 100^ 



Ck = C/100 for 21<k< 100, Ck = for k > 101, 

where C is the normalization factor. In Fig. 4 we plot the error ratios for the wavefunction (A) with 
the second-, the fourth- and the sixth-order decompositions as a function of the effective Trotter 
number N[, which is the Trotter number multiplied by the number of S'2's in Eqs. (p^, ( p^ ) and 
(|l^)P|. Results for the wavefunction (B) with the fourth-order decomposition are also plotted for a 
comparison. We see that the errors from each decomposition decrease as 0{1/N^) in agreement 
with the theoretical expectation. It should be noted, however, that the fourth-order decomposition 
is more recommendable than the sixth-order one because the error from the latter turns out to 
start decreasing only for the large value of Nt. 

Fig. 5 is to demonstrate that the operator e~^^^~^^ can indeed pick up the states with 
the energy ~ E. We calculate L ■ {(f)'^''^ \ e-^'-"-^^ \ (f) {I < k < L) using the wavefunction 
(f>i = + 4>f'^ -I- • • • -f- (j)i"^)/\fL for I 0) and plot them versus Ek — E. We see that the values 
of {(f)^^^ I e~^'^^^^^ I (f)) shape a sharp peak in the narrow range around Ek = E, while they 
are negligible for other fc's. This powerful selectivity should be emphasized to be one of the most 
advantageous features of our method. 

Finally, we show the results on I{E) and P/n ■ I{E) in Fig. 6, where the diamonds are the 
numerical data and the dashed (solid) line is the analytic result with (without) the extra factor ^. 
In the range 200 < (3 < 30000, the Gaussian function multiplied by the extra factor is almost flat 
with the value ^ 0.5, which is the exact value in the continuum limit L 00. This result endorses 
that ^/P/n ■ I{E) fills the role of the delta function. 



6 Three-dimensional problems 

In this section we discuss the problems of three dimensions. Here the Hamiltonian is 

H = -V^ + V{r). (23) 

^Namely Af^' = Nt {N'^ = bNt, = 14A''f) for the second-order (fourth-order, sixth-order) decomposition. 
^Remember that, in order to approximate tlie delta function by the Gaussian function, we need the extra factor 
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In calculating exp[—/3{—W^ + V{f) — E}^] we simply extend the method developed in the one- 
dimensional case and add terms of 

d^ 



dx'^dy'^ ' dy^dz^ ' dz^dx"^ ' 

for which we use the replacements 

d'^4>{x,y,z) 1 

— ^rj.2^y2 -^[(Pi+l,i + l,k + (Pi+l.j-l,k + (Pi-l.j + l.k + <Pi-l,j-l,k 

- '2{4'i+i,j,k + (t'i-ij^k + 4>2,j+i,k + 4'i,j~i,k) + 4:(l}i,j,k], (24) 

and so on. 

Let us consider the three-dimensional harmonic oscillator, whose potential is given by 

V{r} = X\r\^=X{x^+y^ + z^), (25) 

where we also fix A = 1/4 in our calculations. The energy eigenvalues and their eigenstates can be 
obtained from those in the one-dimensional case. The state | 0) is given by the wavefunction 

i,j,k=l 

where represents an eigenstate in the one-dimensional case and we take 



Cimn = l/VlOOO for /<10, m<10, n < 10. 

Fig. 7 plots the numerical results on I{E) = (0 | e^^'-^^^^^ \ 4>) obtained by the second-order and 
the fourth-order decompositions, which indicates that the error is 0(1/ N^) in these cases, too. 
We then apply our method to calculate the density of states defined by 

p{E) ^tr[5{H -E)\. 

Replacing the delta function with the Gaussian function 

we evaluate the trace by means of the random states||l|, 

L 

I = X! {{£,ijk£,i'fk'))=S.H'Sjj>Skk', 

i,j,k—l 

where a state | i, j, k) denotes a particle at the position (x^, yj, Zk) and ((•)) represents the statistical 
average. Results on /3) with our method are shown in Fig. 8 by the dots. The solid line in 
the figure gives the exact value of p{E, 0) calculated using the all eigenvalues in the discretized 
three-dimensional space, 

L 



i,j,k=l 



7 



where E{i, j, k) = Ei + Ej +Ek and E^jj^-^ is the k)th energy eigenvalue for the one-dimensional 
harmonic oscillator. We also show E^ by the dashed line, which is obtained in the continuous space 



where '-Ox{y,z) denotes fca;(j/.z) ^ 1/2- We see that the agreement is satisfactory. 

7 Summary 

In this paper we present detailed descriptions of a new method to calculate expectation values 
of a delta function of the Hamiltonian. We replace the delta function by the Gaussian function 
and apply the Suzuki- Trotter decomposition to it. For concrete examples we carry out numerical 
calculations on the quantum mechanical problems. Our results for the harmonic oscillator problems 
in one- and three-dimensional space indicate that this method is useful to study the dynamical 
quantities. 

Let us close our summary with a few remarks on further possible developments of our method. 
It would be very important to apply it to the problems beyond the quantum mechanics, the spin 
system[H for instance. Applications to the Green function and a step function are also quite 
interesting. 
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Figure 1: A schematic picture for the matrix H2, where crosses (x) indicate its non-zero matrix 
elements. H^^^ {Hf\ Hf \ iff ^) has non-zero elements only within the area surrounded by the 
bold (solid, dotted, dot-dashed) line. 




Figure 2: Errors in our calculations on (0^^) | e~^^^~^^^ \ (p^^^^)) in the fourth-order Suzuki- 
Trotter decomposition with (3 = 100 and E = 0.094716 versus the Trotter number Nt. Here | 0^*'^) 
denotes the fcth eigenstate of the Hamiltonian and the error is defined by the absolute value of 
(^(fc) I e-0("-^f I </,(20)) - e-/'(^=o-B)= . Sk,2o for each k. 
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Figure 3: Errors on {(jP'^ \ e~^^^~^^ \ (f>^'^^^) in the fourth-order decomposition for fc < 40 with 
fixed Trotter numbers Nt = 1000 (diamonds) and Nt = 2000 (squares). Values of the parameters 
(3 and E are the same as those employed in Fig. 2. 
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Figure 4: The error ratio of I{E) with E = 0.094716 defined by (the approximate vahic — the 
exact value) /(the exact value) with the second-, the fourth- and the sixth-order formulae versus 
the effective Trotter number N^.. The value of P is 200 for the wavefunction (A), while it is 2000 
for the wavefunction (B). 
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Figure 5: Absolute values of L ■ (^('^^ | e'^^'^"-^'^^ \ (f)) {I < k < L and E = 0.96302) plotted as a 
function oiEk-E for (3 = 100, Nt = 2000 (the dots) and for /3 = 400 and Nt = 8000 (the pluses), 
using the fourth-order decomposition. 




Figure 6: The numerical data (the diamonds) and the analytic result (the solid line) on I{E). 
The analytic result on (/J/tt)^/^ • I{E) are also presented by the dashed line. Here we use the 
wavefunction (B) with E = 0.094716. The Trotter number Nt is fixed to be 20/3. 
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Figure 7: The error ratio of I{E) versus the effective Trotter number iV/ for the three-dimensional 
harmonic oscillator. Here = 32^, Xmax = —Xmin = 8.0 , /3 = 4.0 and E = 1.1101. 




0.1 1.0 10.0 

E 



Figure 8: p{E, (3) with (3 = 32.0 versus the energy E for the three-dimensional harmonic oscillator. 
Here = 32^, 

Xmax — ~Xmin — 8.0 and Nt — 300. Results of our calculation (dots) are obtained 
by the averages over 10 samples of the random states. 
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